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1. Introduction 

In order to reduce the computational complexity of light beam propagation 
simulation, a number of approximations are typically employed. We take advantage 
of these approximation techniques for the paraxial case to derive an efficient and 
robust method that allows the application of conventional finite difference schemes 
on a uniform grid in the computational space, despite the difficulty of an interface 
present in the domain. 

Prom Maxwell's field equations describing the behavior of monochromatic light, 
we obtain the time-dependent Helmholtz equation. 



where E = E{x,y, z,t) is the electric field intensity, is the Laplacian operator, 
and c is the phase velocity, or speed of light in a particular medium. 
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Let E be the field intensity of a monochromatic plane wave of the form 

E{x,y,z,t) = U{x,y,z)e'^^''', (1.2) 

where u is the frequency of the light. Then from (1.1) we acquire the time- 
independent Helmholtz equation 

{V^ + K^)U{x,y,z)=0, (1.3) 

where k = l-Kvjc = 2-7t/X is referred as the wave number, and \ = c/u is referred 
as the wavelength. Functions \U{x,y,z)\ and aTg{U{x,y, z)) are the amplitude and 
phase of the wave, respectively. Assume that the z is the direction of the beam 
propagation. We may further consider the wave function with a complex amplitude, 
that is, 

U{x,y,z) =u{x,y,z)e-''^^ (1.4) 

where u is called a complex envelope. A paraxial wave becomes realistic if the 
variation of u is slow in the z-direction. 
Substitute (1.4) into (1.3) to yield 

v^u(x, .) - m^^^^^ + = (1.5) 

where 

92 92 
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dx^ dy^ 



V^ = — + 



is the transverse Laplacian operator. In a paraxial case, we may assume that 
within a wavelength of the propagation distance, the change in u is sufficiently 
small compared to \u\ [12]. Thus, 



dz"^ 



which indicates that 

~ 

dz'^ 

Therefore we arrive at an approximation of (1.5), 

V72 / ^ o- 9u{x,y,z) 

Vtu[x, y, z) - 2iK — = 0, (1.6) 

which is called the slowly varying envelope approximation of the Helmholtz equation 
[1, 6, 13]. 
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Under the transformation r = yaP^+y^ and cj) = arctan(y/x), (1.6) can be 
reformulated to 

+ + _2m— ^ u(r z) -0 (17) 

The equation (1.7) has been utihzed frequently in laser beam propagation simula- 
tions in the past decades [6, 7, 16]. 
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Figure 1.1. An illustration of the leans area [15]. 



In this paper, we shall consider a cylindrically symmetric domain for spherical 
lens environments. For the situation we may assume [6] that 







and (1.7) can be simplified to yield 



du 

2zK-(r,.) 



d^u 1 dv> 

-^{r,z) + - — {r,z), < r < ro <C oo. (1.^ 



We assume that the wavelength of light, A, is 9.449 iim and that light is incident 
from air (refractive index ni = 1) into glen (refractive index n2 = 1.5) so that we 
may adopt the following wave numbers [6] 



K(r, z) 



Kq 



27r?ii 2 ■) _i 

" — X 9.97543 X 10 cm , in medium one; 

3 ' (1.9) 

9.97543 X 10 cm , in medium two. 



X 
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The above implies that a coefficient in (1.7) is discontinuous at the lens interface. 
Relation (1.9) represents a single surface situation. Multiple surface scenarios can 
also be discussed with additional k values. Needless to say, the discontinuity adds 
considerable difficulties to the computation of the numerical solution of the differ- 
ential equation [3, 5, 9, 15]. 

We employ Neumann boundary conditions 

Ur{z, 0) = Ur{z, i?i) = 0, Z > 0, (1.10) 

at the bottom, r = and top, r = i?i of the rectangular domain. 

For the initial solution of boundary value problem (1.7), (1.8), we use the 
following approximation of a Gaussian beam with point source [4] 

where Po is the Gaussian beam width, while (5 and A are parameters such that 

i9 = — — = — ^ — A= p'-^^o 



2. Base Difference Scheme and Stability 

Let < z < Z for (1.8)-(1.11). Further, let h = Ri/M, r = Z/N, where M, N e 
are sufficiently large. We may introduce the uniform grid region, 

nh,r = {{mh, nr) \ Q < m < M,Q < n < N} 

over the rectangular domain Q, used. For the sake of simplicity, we denote = mh 
and Zn = nr. In addition, we will use z^-a) < a < 1, for specifying a non-grid 
point between Zn-i and Zn whenever needed. 

Let us start with a linear second-order partial differential equation of the form 

C5^^ + C4^ + C3— -FC2^ -Fciu-l-co = 0, {r,z)&Q,, (2.1) 

together with (1.10), (1.11). The coefficients Cj of (2.1) are functions of z and r 

and may be discontinuous due to (1.9). 

We propose a six-point, two-level Crank-Nicholson type scheme for solving (2.1) 
and (1.9)-(1.11), 

C5 



2/lT 



[^m+l,n f-m— l,n '^m+l,n—l ~l~ ^m— l,n— l] 



"I" [^"^+l>w 2tij7i^n -|- Um—\,n "I" "^m+ljU—l '^'^m,n—l "I" "^^m— l,n— l] 

~^ [^m+l,n Ujn—l,n "I" ^m+l,n— 1 "^^m— l,n— l] 

C2 C'[ 
H ['^m,n ^m,n— l] [^mjn ~l~ ^m,n— l] ^~ Cq — 0. 

7" ^ 
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The above can be conveniently reformulated to a partial difference equation 
\2hT + 2/i2 + Ah) + V /?2 + r + 2 J + I" 



Um—l,n 



2hT 2/i2 4/i. 

/_C5 C£ _ C3 N / C4 C2 _ Ci\ 

l2/ir 2/^2 4/1 J ^'n+l.n-l + ^2 + ^ 2) """^'""^ 
/ C5 C4 C3 \ 

l-2^-2^ + 4^J^"^-^'"-^+^°- ^^-^^ 



Let be a sufficiently smooth function defined on $7. We define 

d^w d'^w dw dw 

Pw{r, z) = C5-Q^ + C4-^(r, z) + c^-^{r, z) + C2-^(r, z) + ciw{r, z) + cq 



and 



/ C5 C4 C3 \ / C5 C4 C3 \ 

+ l"2^ + 2^ " 4^J + l"2)^ + 2^ + "^m+l,n-l 

/ C4 C2 Ci \ / C5 C4 C3 \ 

+ (-^ - - + yj w^,n-i + + ^ - ^-i,n-i + CO. (2.3) 

By expressing the function w at the grid points as Taylor expansions evaluated 
at reference point [rrm ^n-^ substituting these into (2.3), it can be shown 
that 

||(P-P;,,,)u;^,,_i/2|| =0(/i' + /iV + r2). (2.4) 
Let a = r/h be bounded and cr — )■ as /i, r ^ 0. Then (2.4) ensures 

1. the consistency of the finite difference scheme (2.2); 

2. a second-order local truncation error of the scheme (2.2); 

3. the numerical stability depends on the particular coefficients of the differential 
equation considered. 

Definition 2.1. Consider a homogeneous finite difference scheme written as a 
system of linear equations as below: 

BUn = CUn-1 

or 

where vector Un = {uk^n}^=n o^i^d the difference operators B,C E C^^^ are coeffi- 
cient matrices. Let E = B~ C. If there exists a constant > independent ofn, 
h and r such that \\E^\\ < K for some norm \\ ■ \\, we say that the scheme is stable 
in the Lax-Richtmyer sense [10, 11]. 



6 



Gonzalez, Guha, Rogers and Sheng 



Due to the inclusion of a cross-derivative term in our transformed equation, 
our difference scheme wiU not be stable in the Lax-Richtmycr sense. We define a 
notion of practical stability which holds for a range of propagation step sizes that 
afford us sufficient resolution in our simulations. 

Definition 2.2. Let p{B~^C) he the spectral radius of kernel matrix B~^C. If 

p{B-^C) < 1 

at all propagation steps < nr <T with a transverse direction step size eo < h < ei 
for some ei > eo > 0, we say that the scheme 

BUn = CUn-1 

is stable within a parameter range. 

While this stability condition does not specify a norm for convergence, it does 
guarantee that perturbations will not increase exponentially with n. 

Definition 2.3. A matrix A G g^id iq positive semistable if every 

eigenvalue of A has nonnegative real part. 

Tlieorem 2.4. Let A,B,C,G e C^^^ be such that 

B = G + A, C = G-A. 

Then the difference scheme defined by 

BUn = CUn-1 

is stable if and only if G~^A is positive semistable. 

Corollary 2.5. Let A,B,Ce C^^^, d be a positive real number such that 

B = dl + A, C = dl-A. 
Then the difference scheme defined by 

BUn = CUn-l 

is stable if and only if A is positive semistable. 

Recall (1.8). We have the corresponding paraxial Helmholtz coefficients for the 
general equations (2.1) and (2.2), 

1 ^ ^ 

C5 = 0, C4 = 1, C3 = -, C2 = -2iK, C\ = 0, Cq = 0. (2.5) 
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It follows therefore that (2.2), (1.10) and (1.11) can be simplified to the following 
homogeneous paraxial Helmholtz difference scheme, 

-a ^1 + Um+l,n + (2 + 2a) Um,n - a - Um-l^n 

= a ^1 + Um+l,n-l + (2 - 2a) Um,n-1 + a ^1 - Um-l,n-l, (2-6) 

t^m,o = e-'^"^/*, (2.7) 

-2aui^n + (2 + 2a) uo.n = 2q:«i,„_i + (2 - 2a) no,n-i, (2.8) 

(2 + 2a) UM,n - 2auM-i,n = (2 - 2a) UM,n-\ + 2auM~i,n-i, (2.9) 

where 

ri 

a = 



2k/?2 ■ 

Following the analysis method outlined above, we can express our scheme in 
matrix form 

BUn = CUn-1 

where B = G + A, C = G — A, G = 2I, and A is tridiagonal. Investigating 
properties of the eigenvalues of ^ = {am,n}, where 



^ ... 1, M-i, 

2m , 





= 2a, 


Om,m— 1 


-A 


aM,M-l 


= -2a, 




= —a ^ 


ao,i 


= -2a, 



^ ;, m = l,2,...,M-l, 



w are able to show that the eigenvalues of matrix A are purely imaginary, and thus 
have nonnegative real parts. Thus, A is positive semidefinite. By Corollary 2.5, we 
can show the following. 



Theorem 2.6. Let n he a constant. Then the homogeneous paraxial Helmholtz 
difference scheme (2.6) -(2.9) is stable. Further, there is lower boundary restriction 
on the step size parameter, h, in this case. 



3. z-Stretching Domain Transformation 

One possible way of avoiding the computational difficulties presented by the dis- 
continuity of K at the interface is by decomposing the domain into three sections, 
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pre-lens, lens, and post-lens, and stretching each segment by one-to-one transfor- 
mations onto rectangular areas. We would then be able to use conventional finite 
difference techniques, such as that introduced in Section 2, to solve (2.1) together 
with initial-boundary conditions on each segment. The grid stretch can be achieved 
either in the direction of electro-magnetic wave propagation z, or the direction of 
r. Each of the approaches have distinct advantages. We will only focus on the 
former strategy in this paper. In this case, the numerical solution computed at the 
rightmost edge of the pre-lens segment becomes the initial condition of the next 
segment. 

Let r = r{^,Q, z = z{$,X) be the one-to-one stretching transformation to be 
used. Thus, 

du du du d( du du du d( 
dr 8$, dr d( dr ' dz d$, dz d( dz ' 
d'^u dud'^^ ^ d'^u I ^ d'^u djdC, ^ dud'^C ^ d'^u fdC 



dr"^ d^ dr^ d$,^ \dr J d^d(^ dr dr d(^ dr"^ dC,"^ \ dr 
A substitution of the above into (1.8) yields 

2 . 'dud^ ^dudC,\ du d"^^ ^ d'^u V _|_ 2 d^dC 



5^ dz dC, dz J 8^ dr"^ d^"^ \dr J d^d( dr dr 

dud^C d^fdCV i 
d( dr^ d<^^ \ / r \d^ dr d( dr 



which can be regrouped into 



dz dr^ r dr ) dC, \ dz dr"^ r dr J d^, 

\dr ) de \drdrj d^dc \dr J dC ■ ^ ' 

To map the lens area fl = [O < z < Z, (z — -I- < B?, r > 0} into a rect- 
angular area ^ = <C, < Z, 0<^<i?i},a natural choice is the following trans- 
formation, 

^(r, z) = r, C{r, z) = — — -===Z. 



Z - R + VW^ 



In this particular case we have 

^ = 1, ^ = 0, ^ = 0, 

dr dz dr^ 



dC _ rZjz - Z) dC^Z_ d\ _ - Z) [pR^ + 2r'VW^) 
dr pVi?2-r2' dz p' dr^ [R-^ - r'^f'^ 
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where p = Z — R + We define functions 

</'(e,c) ■■ 



R 



^(^,0 := |^(r(e,C),^(C,C)) = 



Z 



2 ' 



Subsequently, (3.1) can be simplified to 



I ,\ du Idu d' 



u 



+ 



u 



(3.2) 



It can be demonstrated that at every point in the transformed lens segment, we 
have 



u 



< 



RjZ 



[R^ - Rf) \z-(r- Vi?2 - r(^ 



d u 



0. 



Thus for a typical lens where the slowly varying envelope approximation is 
applicable, (3.2) can be simplified to a more appropriate form for computations 
within the transformed lens area. 



1 \ du Idu (?2 



u 



dm 



(3.3) 



Lens Segment: Physical Space 




Lens Segment: Computational Space 



Figure 3.1. LEFT: An in-lens domain before a z-stretching. RIGHT: The in-lens 
domain after a z-stretching. 
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In Figure 3.1, we illustrate the z-stretching by showing an example of the grid 
within the lens before and after the stretching. It is evident that the transforma- 
tion enables the six-point scheme to work accurately and efficiently. Instead of 
(2.5), now we have the following set of adjusted coefficients for the in-lens segment 
difference scheme, 

C5 = 2(p, C4 = 1, C3 = ^, C2 = -2iK0 + 1p + ^(/), C\ = 0, Co = 0. 

Subsequently, (2.6) can be modified to 



-7 



26 ( 1 



+ {2 + 2a'y)um,n - 7 



7 



+7 



2(j) ( 1 

T+"i'-2;;^ 



2(/> 



+ a 1 



'^m+l,n— 1 

-I- (2 - 2a-i)Ujn,n-\ 



Itm— l.n— 1 ) 



1 

2m 



'^m— 1,; 



(3.4) 



where 



and 7 = 7(^, C) = 2iKd - - 



In this method, the Gaussian beam input equation is evaluated at the lens 
surface, and becomes the initial solution of the simulation at the edge of the lens 
segment. No computation is necessary in the prc-lcns segment. The solution at the 
right edge of the lens segment becomes the initial solution of the post-lens segment, 
which we simulate using the homogeneous scheme described earlier. 

To determine the boundary conditions applicable within the lens segment, note 
that 



dr 



(^, C) = ^ (^, C) ^ (C, C) + ^ (?, C) ^ (?, C) , 



thus 



^(0,0 = ^W.C) 



0, o<c<z. 



(3.5) 



It is interesting to note that for a lens that tapers to a point at the top, the geometric 
interpretation of this new boundary condition is that the single point (i?i, Z) has 
been stretched into the upper edge of our transformed rectangular domain, i.e. the 
upper boundary in the computational space corresponds to the single point {Ri , Z) 

in the physical space. 

To demonstrate stability, wc use the matrix analysis method introduced in the 
previous section. For our scheme -Bu„ = Cu„_i, we have B = G + A, C = G — A 
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with matrix G = {gm,n} and A = {a^_„}, where 

gm,m = 2, m = 0, 1, . . . , M, 
9m,m-l = -j^, m = 1,2, ... ,M- 1, 

gM,M-i = 0, 

-270 

gm,m+i = ^ , m = 1,2,. . . ,M - 1, 
£^0,1 = 0, 

am,m = 2q;7, m = 0, 1, . . . , M, 

am,m-i = (^-'^ ^ 2m) ' = ••• 1, 

aM,M-i = -2a7, 

am,m+i = (^^ + ' rn = l,2,...,M -1, 

ao,i = -2a7. 

Examining the real part of function 7 (^, ^) resulting from our chosen transfor- 
mation, we see that matrix A is positive semistable. This property will still hold 
if the transformation is adapted to other convex lens shapes. Further, if we have a 
transverse step size h such that 

/i>2|7(e,C)0(e,C)l for 0<^<Ru 0<C<Z 
or equivalently, if the number of grid points in the ^ direction, M, is such that 

where 

hmin = 2 rnax I7 (^, C) (j) (^, C) I 

then matrix G + G* is positive definite. Then based on Theorem (2.4), we can 
prove the following. 

Theorem 3.1. Let k be discontinuous as given by (1.9) and 

/i>2|7(e,C)<^(e,C)l for 0<^<Ri, 0<C<Z. 
Then the difference scheme (3.4) -(3.5) is stable on z-stretched domains. 

For the parameter values utilized in our simulations 

K = 9.97543 X 10^ R = 1.969, Z = 0.7643, Ri = 1.5574 (3.6) 

we need M < 1.2092 x 10^. We list a sampling of maximum grid points M for other 
parameter values. 
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z 


k = 8000 


k = 10000 


k = 12000 


0.1 


7439 


9365 


11279 


0.3 


5642 


7046 


8848 


0.5 


4059 


5064 


6068 


0.7 


2496 


3104 


3711 


0.9 


1033 


1257 


1479 



Table 3.1. Maximum grid points in the transverse direction with R = 1. 
4. Simulation Results and Observations 

All simulated results were implemented on dual processor Dell workstations with 

at least double precision. MatLab, Fortran and C++ programming languages 
were utilized. Dimensionless models are used throughout the computations. For 
the sake of simplicity, we do not tend to re-scale numerical solutions back to their 
original physical dimensions in simulations. 

In the following numerical experiments wc use parameters (3.6) listed in the 
previous section. We further select h = Ri/M where M = 5 x 10^ and r = Z/N 
where N = 1.6 x 10^ in the six-point scheme used in the computational solution 
space after a designated z-stretching. 
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real (r), imaginary (g) and intensity (b) of the initial value 




Figure 4.1. Normalized initial value function u(r, 0). Red curve is for the real part 
and green curve is for the imaginary part of the function. Highly oscillatory features 
of the function is clear. 

We show the real part of the simulated solution, a{z) = real{u(0, z)}, in Figures 
4.2-3. We may observe that while the function value of a is relatively stable before 
the focusing point zj ^ 0.94778, it increases dramatically as z — ^ Zf. This can be 
viewed more precisely in the enlarged picture of Figure 4.3. 




0.1 0.2 0.3 0.4 0.5 



0.7 o.a 



1.1 1.2 1.3 1 .4 1.5 



Figure 4.2. Real part of the simulated solution at the center point r = 0. The 
numerical solution increases rapidly as z approaches the focusing location. Then the 
simulated oscillatory wave diffuses after the focusing point. 
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real part of the solution 




Figure 4.3. More detailed real part of the simulated solution near the focus point, 
r = 0. Same conditions as in Figure 4.2 are used. 

Figures 4.4-5 are devoted to the imaginary part of the numerical solution, b{z) = 
imaginary{ii(0, z)}. Similar to the real part, b is relatively stable before the focusing 
point Zf ~ 2.7431 and is highly oscillatory as z — t- Zf. The phenomenon can be 
viewed more clearly in the enlarged picture of Figure 4.5. 



imaginary part of l 



/ 


i 







O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 



Figure 4.4. Imaginary part of the simulated solution at the center point r = 0. The 
numerical solution increases rapidly as z approaches the focusing location. Then the 
simulated oscillatory wave diffuses after the focusing point. 
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imaginary part of the solution 




Figure 4.5. More detailed imaginary part of the simulated solution near the focus 
point, r = 0. Same conditions as in Figure 4.4 are used. 

Define the numerical intensity function as 

r(r, z) = ^J real^ [u(r, z)] + imag^ [n(r, z)] >0, 0<r<i?i, 0<z<Z. (4.1) 

In Figures 4.6-7 we plot this intensity function against the propagation direction 
z as r being chosen as zero. It is interesting to find that the intensity increases 
rapidly as z — t- zj. The observation is consistent with our previous results. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 OA 



1.1 1.2 1.3 



Figure 4.6. Numerical intensity function of the simulated solution at the center 
point r = 0. The intensity increases rapidly as z approaches the focusing location. 
Then the intensity value diffuses out after the focusing point. 
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Figure 4.7. More detailed numerical intensity function of the simulated solution 
near the focus point (r = 0). Same conditions as in Figure 4.6 are used. 

As a comparison, we further plot the numerical intensity function at the loca- 
tions near and at the focusing point in Figure 4.8. It is observed that the numerical 
estimate of the intensity oscillates rapidly in the r-direction. The intensity increases 
sharply near the center point of the lens, r = 0, while z approaches the focusing 
point location. The simulated wave profiles well match the experimental results. 
The algorithms can be used to provide reference values for further explorations. 

real part of u at focusing position, z = 0. 94778 cm 



0.25 0.5 



Figure 4.8. (ask authors for the image) Highly oscillatory numerical intensity 
function of the simulated solution at the focusing location z = 2.7431. The computed 
intensity increases exponentially near the center, as compared with much lower profiles 
away from the focusing area. 



5. Conclusions 
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In this discussion, we have employed a z-stretching domain transformation to 
map the lens and post-lens domains into convenient rectangular shapes, where 
we can utilize well established finite difference methods on uniform grids. The 
resulting method provides a useful approximation technique that can be efficiently 
implemented with less than 50 lines of Matlab code in the simulation loops. 

A more powerful strategy may be the use of optimally combined z and r stretch- 
ing transformations to increase computational resolution and accuracy of the nu- 
merical solution in critical local regions. Magnifying a particular subregion in the 
transformed coordinate space is computationally equivalent to increasing the refine- 
ment of the grid in that region using techniques such as adaptive mesh refinement 
[2, 14]. Alternatively, domain transformation could facilitate the use of established 
adaptive grid techniques, since applications of standard mesh refinement or moving 
mesh technics are straightforward on rectangular domains [8, 14]. Further, hyper- 
bolic smoothness maps have also been proved to be extremely useful auxiliary tools 
to consider in practical optical beam computations [15]. Detailed discussions will 
be given in our forthcoming reports. 
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